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Abstract. It is well known that increase of the spatial dimensionality enhances the 

fluid-fluid demixing of a binary mixture of hard hyperspheres, i.e. the demixing occurs 

r^ ' for lower mixture size asymmetry as compared to the three-dimensional case. However, 

O , according to simulations, in the latter dimension the fluid-fluid demixing is metastable 

with respect to the fluid-solid transition. According to the results obtained from 

approximations to the equation of state of hard hyperspheres in higher dimensions, 

the fluid-fluid demixing might become stable for high enough dimension. However, 

f~^ ' this conclusion is rather speculative since none of these works have taken into account 

CO . the stability of the crystalline phase (by a minimization of a given density functional, by 

spinodal calculations or by MC simulations). Of course, the lack of results is justified 

by the difficulty of performing density functional calculations or simulations in high 

dimensions and, in particular, for highly asymmetric binary mixtures. In the present 

^) ' work, we will take advantage of a well tested theoretical tool, namely the fundamental 

measure density functional theory for parallel hard hypercubes (in the continuum and 
in the hypercubic lattice). With this, we have calculated the fiuid-fluid and fluid- 
solid spinodals for different spatial dimensions. We have obtained, no matter what 
the dimensionality, the mixture size asymmetry or the polydispersity (included as a 
C^ ' bimodal distribution function centered around the asymmetric edge- lengths), that the 

fluid-fluid critical point is always located above the fluid-solid spinodal. In conclusion, 
these results point to the existence of demixing between at least one solid phase rich in 
large particles and one fluid phase rich in small ones, preempting a fluid-fluid demixing, 
independently of the spatial dimension or the polydispersity. 



PACS numbers: 64.75. Gh, 64.75. Cd, 64.60.De, 64.70.qd 
Keywords: classical phase transitions (theory), phase diagrams (theory), monodisperse 
fluids, mixturs and polydisperse fluids, discrete fluid models 
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1. Introduction 

Entropy-driven demixing was an active line of research since the 1990s. Several 
works had shown the presence of fluid-fluid (F-F) demixing in binary mixtures of 
hard anisotropic particles [H El |3l HI |5l E]. The demixed phases can have different 
symmetries, with or without long range orientational order, the so called nematic or 
isotropic phases respectively. However, the simplest binary mixture of hard particles 
one can imagine, namely the hard sphere (HS) mixture, does not demix into two fluid 
phases even for high mixture size asymmetry. The first simulation study on the fluid- 
solid phase separation in the binary mixture of hard spheres |i7] does not report any 
F-F demixing. However, they only consider mixtures with diameter ratios very close 
to one. More recent theoretical and simulation works on this system consider more 
asymmetric mixtures and for large enough asymmetry they have shown that the F-F 
demixing is metastable with respect to the fluid-solid (F-S) phase separation, with the 
crystalline phase rich in large spheres while the fluid phase is mostly populated by small 
spheres [H [9]. In these works the F-F instability is taken into account in the computer 
simulations by approximating the effective hamiltonian through an explicit depletion 
potential approximation (integrating out the degrees of freedom of the small spheres). 
The same scenario occurs in binary mixtures of parallel hard cubes (PHC) with particles 
freely moving over space ^U\ or constrained on a lattice [H]. Note that although the 
cubes are anisotropic particles, the gain in the entropic volume after demixing can 
only be reached by changing the mixture composition or the density. As the cubes are 
parallel, the orientational entropy does not play any role here. In this sense the mixture 
of PHC is similar to the HS mixture. However, the effective depletion pair potential 
between two large solute particles immersed in a solvent of small particles is greater if 
the particles have cubic symmetry. As a matter of fact, this is why the scaled particle 
theory (SPT) applied to cubes predicts a metastable F-F demixing ^0\, while the same 
theory applied to HS does not. 

In contrast to mixtures of additive HS, mixtures of non-additive HS can exhibit 
a stable F-F demixing as was shown by several authors both theoretically and by 
simulations [12]. This behavior can be perfectly understood by resorting to the 
following argument: when the interactions between different species are more repulsive 
(with a larger excluded volume) than those corresponding to the same species, the 
particles in the demixed configurations have a lower entropic volume and thus a greater 
configurational entropy (and consequently a lower free-energy). 

The study of the hard hypersphere (HHS) one-component fluid has received great 
attention in recent years because certain questions that are ambiguous in 3D, e.g., 
jamming, crystallization, glass formation, find clearer answers when examined from a 
dimensional perspective. Also the poor knowledge of the lattice packing properties of 
HHS in high dimensions has led to works toward that direction [13]. Recent works 
have shown that the freezing of HHS is frustrated because configurations of spheres 
with simplex liquid ordering are very different from those corresponding to the periodic 
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lattice geometry [14]. Also calculations of the pair correlation functions g2{r) evince that 
the short-range ordering decreases appreciably with dimensions \T5\. The frustration 
translates into higher free-energy barriers between the fluid and crystal states [HI [15] 
and thus the fluid has a greater propensity (than in three dimensions) to form a 
glass upon compression [15]. Both, simulations [HI [151 [IB] and density functional 
calculations [17] have shown that the F-S transition occurs at lower packing fractions as 
the dimensionality increases and also that the character of this transition remains of first 
order. This does not mean that freezing is favoured when the dimensionality increases, 
since close packing is also reached at lower densities. Finally, the HHS fluid has been used 
as a test for different theoretical approximations of the structural properties of the fluid, 
comparing the direct correlation function c{r) and the radial distribution function g{r) 
obtained from these theories with those obtained by MC and MD computer simulations 
[T81[I91[20]. 

Returning to the demixing behavior of binary mixtures composed of particles 
interacting with hard-core potentials, we should say that the number of theoretical 
studies carried out in dimensions higher than three are very small, to our knowledge there 
is only one [21] . In this work the authors have shown that the F-F demixing of a binary 
mixture of HHS is favored as the dimensionality increases, i.e. the demixing shows up at 
lower mixture size asymmetry when the dimensionality increases. Notwithstanding, this 
is not conclusive enough to discard that one of the demixed phases could be unstable 
with respect to the freezing transition. The only way to settle this question is to involve 
the crystalline phase of a binary mixture in the theoretical or simulation studies, which 
obviously becomes a much more difficult task as compared to three dimensions. 

The aim of the present work is to get some insights into what is the most likely 
scenario: F-F versus F-S demixing in a mixture of parallel hard hypercubes (PHHC). 
We have selected this system because the freezing of PHHC is much less affected by 
the geometric frustration present in HHS and thus the study is easier to carry out. 
For example, the simple hypercubic lattice is probably the most stable configuration of 
particles in a crystalline phase, and the same should occur in the crystalline phase of 
a binary mixture when small cubes occupy a very small fraction of the total volume. 
The other reason behind this selection is related to the theoretical tool we will apply to 
this study: the fundamental measure density functional theory developed for a mixture 
of PHHC defined on a lattice [22] or on the continuum [23]. As has already shown, 
the fundamental measure theory (FMT) applied to a HS fluid gives the most accurate 
results for the freezing transition as compared to other density functional theories [23]. 
Comparison between theory and simulations also confirms the high performance of this 
theory in predicting the equation of state of the crystalline phase of PHC [TO]. Also, 
we have already shown that in three dimensions the F-F demixing is metastable with 
respect to the F-S demixing even for high mixture size asymmetry [10]. Simulations 
carried out on binary mixtures of parallel hard squares found strong clustering of large 
squares with crystalline structure in coexistence with a fiuid of small squares [25] which 
would confirm this scenario in small dimensions. Note that this is in accordance with 
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SPT prediction for this system of a stable mixed phase with respect to the F-F demixing. 

In the present paper we use FMT for PHHC to perform a F-S bifurcation analysis 
which together with the F-F spinodal calculations allow us to conclude that the F-F 
demixing is preempted by the freezing of at least one of the segregated phases (that is 
rich in large hypercubes). We show that this result holds up to the 9th dimension. For 
dimensions greater than 9 the third virial approximation is accurate enough and we can 
use it to show that the same scenario remains up to dimension 25. For higher dimensions, 
we have carried out an asymptotic analysis of the second virial approach (which becomes 
exact in the limit of infinite dimension) and the conclusions do not change: the F-S 
demixing preempts the F-F demixing. Finally, we have included polydispersity in the 
edge length of the PHHC by means of a bimodal distribution function and again we 
find the same scenario. 

The paper is organized as follows. In Sec. [2] we present the expressions for the 
direct correlation functions of a mixture of PHHC as obtained from the FMT defined on 
a lattice and on the continuum. These functions allow us to perform an F-S bifurcation 
analysis and to calculate the F-F demixing spinodals. In Sec. |3] we show the results 
from these calculations. This section is divided into three parts: in Sec. 13.11 we present 
the results from the F-S bifurcation analysis of the one-component fluid of PHHC as 
a function of dimensionality; Sec. 13.21 is devoted to a detailed study of the demixing 
transitions as a function of mixture size asymmetry and dimensionality; and finally in 
Sec. 13.31 we present the results from the bifurcation analysis of a bimodal polydisperse 
mixture of PHHC. In Sec. Hlwe carry out an asymptotic analysis (for large dimensions) 
of the second virial approximation which allows us to conclude that at infinite dimension 
the F-F demixing is metastable with respect to the F-S demixing. The conclusions are 
drawn in Sec. [5] and theoretical details of the bifurcation analysis are relegated to an 
appendix. 

2. Model 

The model consists of a binary mixture of PHHC embedded in dimension n with species 
i = 1,2 having an edge length equal to ai. We define mixture size asymmetry by the 
coefficient k, = cr2/o"i > 1. The mixture is fully characterized by the density profiles 
Pi{r) where the vector r = {xi,X2, ■ ■ ■ , Xn) denotes the position of the center of mass of 
a particle in the n-dimensional space. If the particle positions are constrained to be on 
a simple hypercubic lattice the variables Xi are set equal to i G Z (i.e. we are using the 
lattice mesh size as the unit of length). Note that when the lattice mesh size goes to 
zero, the lattice model goes to the continuum model, therefore in the lattice model we 
have chosen (J2 = 2 in order to emphasize its discrete nature. We will use the density 
functional (DF) based on the FMT obtained by us already on a lattice [22] and on the 
continuum [23]. We refer the reader to those references for a detailed description of 
the way to obtain the DF from the free-energy density of a so called zero-dimensional 
cavity, which can accommodate one particle at the most. For a more general reference 
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on the theory see Refs. 

From the DF we can obtain the Fourier transform of the direct correlation 
functions which are the only expressions we need to compute the F-F and F-S spinodals 
corresponding to the F-F or F-S demixing transitions. From these, we can analyze the 
relative stability between F-F and F-S transitions. 

2.1. Direct correlation functions following the FMT for PHHC 



As we have already shown in [221 [23] the excess part of the free-energy DF corresponding 
to a mixture of PHHC in dimension n, pj-'ex[{pi}], can be obtained from the excess free 
energy of a zero-dimensional cavity 

rj being the packing fraction of the cavity. 

In particular, for the continuum c-component mixture of parallel hard 
hyperparallelepipeds we have 
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ri{r) being the local packing fraction of a maximal zero-dimensional cavity centered at r 
and a^ being the edge length of component i in the / spatial direction. The analogous 
expressions for the discrete case are 

<|..(r) ^ (nAj Mv^'-Hr)), A,f{k) ^ /(I) -/(O), (4) 
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1]^^' (r) being the packing fraction of a zero-dimensional cavity around r. For a connection 
between the difference operator in equation (jlj) and the Mobius inversion formula, the 
reader is referred to [26] . 

Using Eqs. ([2]) and (jlj), we can calculate the direct correlation functions defined 
through the second functional derivative, q" (ri2) = — (5^/3 J-'g^"7(5pj(ri) 5pj(r2), evaluated 
at uniform density profiles Pi{r) = pi. If we calculate the Fourier transforms of these 
functions, evaluate them at the wave vector q = g(l, 0, . . . , 0) and take cr| = ai for any 
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k, in the continuum we obtain 

n-l 



k=0 



c:.^(g) = Ec'r'M.)'^r'" 



X [Ck(Tijjiiq(Tij/2) + Ck+iCTiCTjji{qai/2)ji{qaj/2)] , (6) 

with aij = ai + aj, CJ^ = n\/k\/{n — k)\ the usual combinatorial coefficients and ji(a;) the 
spherical Bessel function of ffist order. The functions Ck{^n-ki ■ ■ ■ ■, ^n) with k = 0, . . . ,n 
can be obtained recurrently from Co(^n) = 1/(1 ~ ^n) using the following relations 

Ck+l{^n-k-l,- ■ ■ ,^n) = 2_^^n-i-l -^ '- , fc = 0,...,n — 1 (7) 



j=0 
.k 



din- 



where we have defined ^k = ^iPi'^i- Note that the packing fraction rj is just equal to 
C,n- For the lattice version we have 

n— 1 ^ 

^ sm (g/2) 



sin(gcrj/2) sin(gcrj/2) sm{q{ai — l)/2) sin(g(crj — l)/2) 



(8) 



1 - ^7/0+1 1 - % 

where we have defined r]k = X]j Pi^ii^^i ~ 1)""'^ for k = 0,1, . . . ,n, which is the uniform 
version of ^^'^''•'(r) with k = ki -\- ■ ■ ■ -\- kn- 

The fluid pressure of the mixture in the continuum is related to the function 
Cniio, ■ ■ ■ ,in) through the relation (n = dPp/d^n- The lattice version give us the 
following expression for the pressure 

n 

Pp = -Y.^i{-ir-'H^-vk). (9) 

fc=0 

We should note that the expressions (E]) and (IHl) in the low density limit pi — )■ 
recover the exact second and third virial approximations. 

2.2. F-F and F-S spinodals 

The packing fraction rj (note that r/ = ,^„ in the continuum and ?7 = ?]„ in the lattice) 
and the wave number q at the F-S transition of the binary mixture can be calculated 
from the singularity of the matrix with coefficients defined through the partial structure 
factors 

Sij {q, v) = 5ij - p^/x'x'jcf^'^ (g) , (10) 

with Xi = Pi/ p the mixture composition and p = J2i Pi the total density. The singularity 
condition gives us, at a fixed composition xi {x2 = 1 — Xi), 

1 - pxicSi^(g) -py/xiX2c']^\q) 
-py/xiX2c(!^{q) 1 - pX2cf^{q) 



A(g,r/,xi) = det , ^ ^^(„) ^(„) 



1 - Pic^ihq)) (l - P2cSHg)) - p,p2 (cS^^g))' = 0. (11) 
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Now, we should look for the lowest value of t] that satisfied this equation. If we take into 
account that A{ri,q,Xi) = defines a function rising{q,Xi), in order to find the values 
^bif(a^i) and gbif at the F-S transition we need to solve 

^bif(a;i) = min rising (g,a;i), 

'" , , (12) 

gbif = argmin?7sing(g,xi). 

ij>0 

In the continuum this is equivalent to the set of equations 

dA{q,r],Xi) 
A(g, r], xi) = = 0, (13) 

In the lattice version, due to the discrete nature of q (in units of the lattice mesh size) 
we do not have the second equality in ( IT3|) . Finally, in both cases the F-F spinodal can 
be calculated from the condition A{0,ri,xi) = 0. 

3. Results 

This section is divided into three parts. The first one is devoted to showing the 
performance of our theory concerning the calculation of the F-S transition of a one- 
component fluid of PHHC as a function of dimensionality n. In the second part we 
analyze the relative stability between the F-F and F-S demixing, again as a function 
of n, of (i) a binary mixture of PHHC and of (ii) a bimodal polydisperse fluid. For 
the former system, we have carried out the analysis both in the continuum and in the 
lattice cases, while for the latter only the continuum system is studied. The last part 
of the section is devoted to the study of the binary mixture phase behavior at infinite 
dimension, by means of an asymptotic analysis. 

3.1. One component fluid 

As we have already shown using both versions of the FMT, the continuum parallel hard 
cube fluid in dimension three exhibits a second order F-S transition where the cubes 
crystallize in a simple cubic lattice [23]; in the lattice version [22] the system shows a 
second order transition where the structure of the ordered phase depends on the size of 
the cubes (from smectic for the smallest cubes possible to solid for larger sizes). Herein 
we will be referring to F-S demixing for both the continuum and lattice model, but we 
ask the reader to take into account that in the discrete case the structure of the ordered 
phase could be different from a solid. In the continuum the same behavior was also 
predicted from simulations [281 EHl [30] . Although the character (second order vs. first 
order), and the symmetry of the crystalline phase were adequately predicted from FMT 
as compared to MC simulation results, the precise location of the transition point is 
underestimated by the theory [10]. This discrepancy is explained by noting that the 
uniform limit of the FMT recovers the SPT and, therefore, the fluid phase predicted 
by FMT inherits all its defects. While all the virial coefficients predicted by SPT are 
positive, it is known that this is not the real case [311 132] where the presence of negative 
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Figure 1. Packing fraction ry (a) and period of the crystalline phase d/a (b) at 
bifurcation as a function of the dimensionality n. 



coefficients results in a poor convergence of the virial expansions. As a consequence, 
although FMT recovers the low-density limit, the equation of state for the fluid near 
the transition point is not very accurate. Notwithstanding, because of the dimensional 
crossover property of the FMT DFs, the predictive power of FMT is greatly enhanced 
when the packing fraction of the crystalline phase is high enough. Actually, for packing 
fractions rj > 0.5 the comparison between theory and simulations is fairly good [TO] . 
Finally, for dimensions higher than three we should remark that there is a rigorous 
result that shows that at infinite dimension the PHHC fluid exhibits a second order 
phase transition to a crystalline phase with a simple hypercubic symmetry |33j . 

We calculate here the F-S transition of the continuum PHHC as a function of the 
dimensionality using a bifurcation analysis, i.e. by solving the equations ( iT3l) in the 
one-component limit to find the packing fraction ?7bif and the period in units of the edge 
length d/a at the transition point (note that d is related to gbif as d/a = 2n/ qui)- These 
values are plotted in figures [1] (a) and (b), respectively. As is already known, the second 
virial approximation becomes exact in the limit n — )■ oo. Here, the packing fraction 
at the transition scales with n as r^ ~ 2^". With the aim of comparison, we have also 
plotted in figure [1] the results from this approach. As we can see from the figure, the 
packing fraction r^bif and the period d/a obtained from both FMT and second virial 
approximations are very similar when n ~ 10. Therefore, it should be expected that 
around this value of the dimension and for higher ones FMT predictions should not 
differ too much from those obtained by means of a virial expansion. This supports our 
approach in the next subsections where we have analyzed the binary mixture stability 
using FMT only up to n = 9 and from there to n = 25 we have used the third virial 
expansion instead, and we have concluded with an asymptotic analysis of the second 
virial expansion, which should be accurate enough for higher dimensions. 
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Figure 2. (a) F-F (dotted) and F-S (dashed) spinodals in the paclcing fraction- 
composition plane for a binary mixture of PHHC in dimension four and mixture 
asymmetries equal to k — 5 and n — 10. The coexistence F-F binodals are also 
plotted with solid lines. The solid circles indicate the position of the critical points 
while the diamonds joined by straight lines indicate some examples of coexisting fluid 
phases (note that one or both of them are metastable with respect to the freezing 
transition). The solid squares indicate the F-S transition of the pure one-component 
fluid, (b) The Gibbs free energy per particle in reduced units minus a straight line 
(to make the coexistence more visible) g* = PG/N — 7.4962 — 1273. 2a;i as a function 
of the molar fraction of the small component. The mixture size asymmetry is k = 10 
and the pressure is fixed to I3p<j\ = 939.49 [coinciding with that corresponding to the 
coexistence points labeled with filled diamonds in (a)]. The filled diamonds indicate 
the F-F coexisting points while the open circle indicates the F-S bifurcation point. To 
the left of this point the fluid branch (with the dashed line) is unstable with respect 
to the solid branch sketched (not calculated) with the dotted line. 



3.2. F-F versus F-S demixing in binary mixtures of PHHC 

In this section, we will analyze the phase behavior of the binary mixture of PHHC as 
a function of the dimensionality. As was already shown in dimension three, the F-F 
demixing is always unstable with respect to the F-S demixing in which a crystalline 
phase preferentially populated by large cubes coexists with a fluid phase rich in small 
cubes [lO]. The uniform mixture becomes unstable with respect to the F-F demixing 
for a mixture size asymmetry k = h + a/24 f« 10 in the continuum and /t = 13 in the 
lattice. However, the position of the critical point of the F-F spinodal curve in the 
packing fraction-composition plane is always located above the F-S spinodal curve, the 
latter computed from the set of equations (fT2l) . This is equivalent to saying that at least 
one of the coexisting uniform phases is unstable with respect to the freezing transition. 
We have found the same scenario for higher dimensions. The main differences are 
only quantitative and can be summarized as follows: the lowest mixture size asymmetry 
at which the F-F phase separation shows up (see table [1]) decreases dramatically with 
the dimensionality and the same happens with the packing-fraction value at the critical 
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n 


3 


4 


5 


6 


7 


^cont 


9.90 


3.73 


2.54 


2.09 


1.87 


^latt 


13 


5 


3 


3 


2 



Table 1. Critical aspect ratio k* for which the F-F demixing appears for the first tirae 
in the continuum and lattice models. 



point. Thus the demixing transition is enhanced with n. On the other hand, the packing- 
fraction values at the F-S spinodal also decrease with n, and although they get closer 
to the F-F spinodal in the neighborhood of the critical point, they remain always below 
it. This behavior is present even for high values of the mixture size asymmetry. It is 
worth noting that n* for the continuum model decreases as the inverse of the dimension; 
actually from the values shown in tabled] we have found that k* follows very accurately 
the law K* = 1 -\- a/{n — b) where a = 3.809 and b = 2.566. The conclusion is that, no 
matter what the dimension or the mixture size asymmetry, the most probable scenario 
is the presence of demixing between a fluid phase rich in small cubes and a crystalline 
phase mostly populated by large ones. This scenario is illustrated in figure [2] where we 
have plotted both spinodals (F-F and F-S) and the F-F coexistence binodals obtained 
from the continuous model of PHHC in dimension four and asymmetries fixed to k = 5 
and K = 10. Note that we have used to represent the composition of the mixture the 
volume composition xi, which is related to the molar fraction as Xj = cr"xj/(cr"), where 
(o""") = J2i^i^i- -^^ ^^ '^^^ ^^^ from the figure, the critical point in the F-F spinodal 
(55, rj'^) is located above the F-S spinodal, therefore at least one of coexisting fluid phases 
is unstable with respect to the freezing transition. This scenario is illustrated in figure 
Mjo), where we plot the Gibbs free-energy per particle as a function of the molar fraction 
of small cubes xi for a fixed value of the fluid pressure (note the difference between xi 
and Xi due to the huge difference in volume between the large and small cubes). From 
the left of the F-S bifurcation point (open circle), the fluid branch is usually located 
above the solid branch which is sketched with a non-convex curve. From our earlier 
calculations in n = 3 the convexity of the solid branch is usually recovered when the 
volume fraction of big cubes is near unity. As we have already shown in [M] the infinite 
mixture size asymmetry limit corresponds to a one-component fluid of sticky cubes. For 
this case, at any fixed solvent fugaticy, the coexistence is between an infinite dilute vapor 
phase and a close-packed crystal. When we include a small amount of polydispersity 
in the edge length of the large cubes the convexity of the solid branch is recovered at 
high packing fractions, and the coexistence is now between a highly dense crystal and 
a vapor phase. We expect a similar behavior for a highly asymmetric mixture. 

In figure El we show the phase behavior of both continuum (figures (a) and (b)) and 
lattice (figures (c) and (d)) PHHC for n = A and n = 5 for a large set of mixture size 
asymmetries, some of them being as high as k = 100. It is remarkable that even for very 
high values of k the behavior is the same as that described above: the F-S demixing is 
the most likely scenario. 
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Figure 3. Phase behavior of a binary mixture of PHHC in the continuum [(a) and 
(b)] and in the lattice [(c) and (d)] for dimensions 4 [(a) and (c)] and 5 [(b) and (d)]. 
The solid and dashed lines represent the F-F and F-S spinodals for different values 
of the mixture asymmetry k (correspondingly labeled in the figure) . The filled circles 
indicate the locations of the F-F critical points. 



This trend can be better appreciated in figures |l](a) and (b) (for tlie continuous 
model) and (c) and (d) (for tfie lattice model), where we plot the coefficient Rn = 
iv^ — ''7bif(a;i))/''7bif(a;i) which measures the relative distance between both spinodals at 
the F-F critical point as a function of n. As we can see from the figure, this coefficient 
decreases with n up to a certain dimension (the position of the minimum) which depends 
on K, and further increases with n diverging at ra — )■ oo (this will be shown by means of 
the asymptotic analysis presented in sectionH]). In these figures are also shown the period 
d/a2 of the crystalline mixed phase with xi = xl calculated at bifurcation as a function 
of n, and the period of the one-component crystal. In general, for low dimensions the 
crystalline phase of large cubes becomes more packed by adding an small amount of 
small cubes, an effect clearly related with the entropic depletion forces. 

Using the third virial approximation [equations (jH]) and ([H]) up to first order in 
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Figure 4. The coefficient i?„ (solid line) and the lattice periods d/a2 corresponding 
to the crystalline mixed phase (the circles joined with dashed line) and to the one- 
component crystal (squares joined with dashed line) as a function of n as obtained 
from the continuum [(a) and (b)] and the lattice [(c) and (d)] models. The values for 
the mixture asymmetries are k — 10 [(a) and (c)] and k = 20 [(b) and (d)]. 



density] we have extended the spinodal calculations up to dimension 25. In figure ([5]) 
we plot the packing fractions of the F-F critical point and of the F-S transition at x^ 
as a function of n in both the continuum and the lattice models. Again, we find the 
expected behavior: the freezing preempts the F-F demixing. 



3.3. Effect of polydispersity on the F-F and F-S demixing 

We have also studied the effect that the polydispersity has on the demixing behavior 
of PHHC We have relegated to the appendix the formalism that leads to the main 
equations for the F-F and F-S spinodals. If we introduce polydispersity in the cube- 
edge lengths around their mean values Ui = 1 and (72 = 50 (i.e. following a bimodal 
distribution function) for n = 3 and (72 = 20 for n = 4, we find that the crystalline phase 
destabilizes with respect to the fluid phase as the polydispersity increases. However, 
although the F-F critical point gets closer to the F-S spinodal it is still located above 
it, even for large polydispersities (see figures [6] (a) and (b)). 
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Figure 5. The packing fractions corresponding to the the F-F critical point and to 
the the F-S transition at x\ as a function of dimensionality following the third virial 
approximation as obtained from the continuum (a) and lattice (b) models. 
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Figure 6. (a) F-F and F-S spinodals of PHHC in dimension three around the F- 
F critical point for different values of the polydisperse coefficient r as labeled in the 
figure, (b) The whole spinodals in dimension four at polydispersity r = 0.577. The 
mixture size asymmetry is fixed to k = 50 in (a) and k = 20 in (b). The circles specify 
the position of the critical points. 



4. Asymptotic analysis 

In this section we will carry out the asymptotic analysis for the functions 7]'^{n) and x\{n) 
corresponding to the position of the F-F critical point as a function of dimensionality 
when n — > oo. For this, we will use the second virial approach, as it becomes exact at 
this limit. The aim of the analysis is to find the relative position of the F-F critical 
point with respect to the F-S transition at large dimensions. In order to do that, we 
will use the asymptotic expressions for the F-F critical point to evaluate ^{qm_in,'n^,x\) 
(with gmin the position of the absolute minimum of A(g, rj'^,x\) as a function of q). Note 
that if this quantity is negative, this means that the F-S transition is below the F-F 
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critical point. For the sake of simplicity, we will restrict this analysis to the continuum 
model. In the lattice case, the expressions become too cumbersome, but qualitatively 
there is no difference with the continuum model, so we have preferred not to include it 
in this paper. 

The expression for the direct correlation functions Cj^ (q) in the second virial 
approximation can be obtained from equation (|6]) as a Taylor expansion of zeroth order 
in the densities. In this approximation, the pressure can be cast as 



t3p = P+T^^PiPj^^^j- 



iU) 



*J 



In the following, we will use the volume composition Xi and the scaled packing 
fraction 770 = ?72" as the adequate variables to evaluate /3p and A(g, 77, xi) when n ^ 1. 
In these new variables their expressions become 

p*{r]o,xi) = Pp{2air 



^0(^1 + - +y 



Z.2 



Xi H h 2a;ia;2 



(15) 



for the pressure and 

A(g,?7o,Xi) = 1 + Xir]ojiiqcri) + X2r]ojiiq(T2] 



Ji{q(^i)ji{q(^2] 



K + K, 



xiX2r]Q 
-2 



4 



Ji(g^i2/2) 



Note that X2 = 1 — xi. The above expression evaluated at g = gives 
A(0, r^o, xi) = 1 + r^o - xiX2Snr]o, 



(16) 



where we have defined the coefficient 6„ 
A{0,rio,xi) = defines the F-F spinodal: 



K-%K + l)/2] 



2ra 



1 + a/1 + 46nXiX2 



1. The equation 



:i7) 



2XiX2Sn 

The critical point can be found as the value of xi at which dp*{riQ^, xi)/dxi = 0. When 
n — 7- 00, the composition at the critical point x^ goes to zero (and therefore X2 x 1); 
taking this into account and using the new variables s„ = \/l + 4(5„xi — 1 and r„ = (5„/k"' 
we have 

dp* 2 + Sr. 



dxi 2s^(l + s^ 



[s^ - 4^/^(1 + y^)s„ - 8r„] , (n^oo). (18) 



The asymptotic solution of dp*{ri^ ,Xi)/dxi = with respect to s„ gives us s^ 
\/t^ + \/ 4\fTn + Tn and then, at the F-F critical point. 



x\{n) X (25, 
and from equation ( TT7|) 



'i\fTn^rn + (1 + v^)\/4v^+^> 



rilin) 



^1/4 



1 + 



1 
2' 



m — 7- 00 



(n^cx)), (19) 



(20) 



From here we have that the leading terms in the asymptotic behavior at the F-F critical 
point are x\{n) x \2^2^j{^K + 1)^/^] and T(^{n) x [2/ (1 -|- k"-^)] , respectively. Note 
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that although ?7o — > oo when n — > oo the packing fraction 77 = 2~^riQ goes to zero when 

n — 7- 00. 

The function A(g, ?7q,x5) is asymptotically 

A(g,r7S,55) >< l-j!iqau/2) + Vo [ji(^^2) - J?(g^i2/2)] , (n ^ 00). (21) 

When n — )■ 00, r^Q — )• 00, and the dominant contribution to A{q,riQ,x1) is just the 
term proportional to ?7§ in equation ( 1211) which for k ^ 1 has its absolute minimum at 
(lmm(^2 = 3.506 with the value A{q^i^, tjq, x^) = — 0.416?7q. This negative value implies 
that the mixture with composition and packing fraction corresponding to the F-F critical 
point is unstable with respect to the F-S transition. Thus, the F-F demixing is unstable 
with respect to the F-S demixing. 

5. Conclusions 

In the present work we have made the first attempt to answer to the question about 
how the dimensionality and mixture size asymmetry affect the relative stability of the 
F-F demixing of PHHC against the F-S phase separation, based on the bifurcation 
analysis for the latter transition. Although the results presented here strongly suggest 
the F-S demixing as the most likely scenario independently of the dimensionality and 
of the mixture size asymmetry, the final confirmation of this phase behavior can be 
given only after carrying out MC simulations or DFT coexistence calculations on this 
system. Both tasks become much more difficult as the mixture size asymmetry and 
dimensionality increase. The main difficulty in performing MC simulations consists 
of the high probability of overlapping many small cubes when a given large cube is 
moved by a MC step. Thus a more sophisticated simulation scheme like a cluster 
movement algorithm is required to carry out these simulations. Even using these 
schemes, the increase of spatial dimensionality (apart from the high mixture asymmetry) 
constitutes an additional complication due to the huge increase in the number of particles 
necessary to perform a simulation with a reliable statistics. Instead of performing an MC 
simulation on the binary mixture, it is possible to map the mixture onto an effective one- 
component fluid with particles interacting via an effective depletion potential which can 
be previously calculated using the formalism described in [8]. Although this mapping 
is approximate it should work very well for very asymmetric mixtures. 

The extrapolation of this result to the HHS mixture or other additive hard mixtures 
should be made with some caution. The elevated barrier for freezing in higher dimensions 
might allow for the metastable F-F separation to be observed and particularly long lived. 
Also, as we have pointed out in section [H the freezing transition in HHS becomes more 
frustrated with dimensionality. Thus it might possibly occur that at some dimension 
this frustration makes the F-F demixing the stable one. To elucidate this question we 
should wait for future works which we hope will have been motivated by the present 
one. 
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Appendix A. Polydisperse PHHC 

In this appendix, we present the formalism used to prove that the transition to the 
crystalline phase preempts the F-F demixing even in polydisperse mixtures of PHHC. 
For this purpose we select a bimodal length distribution function as 

h{a) = xho{a/ai) + (1 - x)ho{a/a2), (A.l) 

where < a; < 1 is the polydisperse composition variable, and we use for /io(o"/o'j) the 
normalized to unity Schultz distribution function with mean values fixed to (Tj {i = 1, 2): 

hoia/a,) = ^^^p^^ f-V exp [-{u + l)a/a,] , (A.2) 

ad (z/ + Ij \aij 

where T{x) is the Gamma function and u controls the width of both peaks located in 

the neighborhood of (Jj and related with the polydisperse coefficient r of the distribution 

function ho{a) as r = \/{o'^)hoia)/{o')h (a) ~ ^ ~ VV^ + 1- The moments of h{a) can 

be calculated as 

Wlkia) = 1^°^ daa-h{a) = Jl\^^^^l^^ -^^ (^< + (1 " ^M) ■ (A.3) 

Near and above the F-S bifurcation point the density profile of a polydisperse 
mixture can be approximated as p{r,a) ^ Po(<^) + f(i"?<^)? where po(<^) = Poh{a) 
is the density distribution corresponding to the uniform parent phase (with po its 
number density) at bifurcation and e(r, a) is a small non-uniform perturbation. The 
minimization of the grand potential 

^[p] = -^id[p] +-^cx[p] - dr dafio{a)p{r,a), (A.4) 

/3Mp] = / ^^ / ^^Z^^^' ^) {ln(V(a)p(r, a)) - 1} , (A.5) 

(where Jid[p] is the ideal part of the DF while po('^) is the fixed chemical potential of 
species with length a and V{a) is the thermal volume of species with length a) with 
respect to p(r, a) gives us 

p(r, a) = po{a) exp {-^^^ + ^f^Ti^) } , (A.6) 

/3ii^{a) = f3fio{a) — ln(V(o')po(o')) being the excess chemical potential of species with 
length a. Near the bifurcation point, expanding the exponential around po(o') we obtain 
an integral equation for the perturbation e(r, a) 

e(r, a) - po(a) f dr' I da'c^""^ (r, a, a'W , a') = 0, (A.7) 



Fluid-fluid versus fluid-solid demixing in mixtures of parallel hard hypercubes 17 

where c^'^\r,a,a') is the direct correlation function of the polydisperse mixture. The 
Fourier transform of Eq. flA.7|l for the wave vector q = g(l,...,0) give us 

e{q,a)-poh{a) f da'c^^\q,a,a')e{q,a') = 0. (A.8) 



The function c'^"^(g, a, a') can be obtained from ([6]) by replacing ctj and aj by a and a' 
respectively. We define the functions 

ef\q) = f daa'uk{q<j/2)e{q, a), A; = 0, 1, (A.9) 

where we have used the notation uo{qa/2) = cos{qa/2), and ui{qa/2) = 2sm{qa/2)/q. 
Inserting c^"\q,a,a') from equation (Q into equation (]A.8|) . integrating over a', 
multiplying the result by a^Uk{qo'/2) {k = 0, 1) and integrating again over a we obtain 
the following linear system with respect to the functions e\ (g): 

n-l 

e?\<l) + Po E (7;f ^^g)er (g) + 7;f ^^n^)ep)(g)) = 0, (A.IO) 

i=0 

where we have defined the matrix coefficients 

7;f -^g) = c^r^ EqomiY^-.-i(g), (A.ii) 

j=0 

r!t''\q) = cr'Y^c] [cMnS^i-U<i) + c,^i^^nili-U<i)\ ■ (A.12) 

j=0 

and the generalized moments 

mf'"^(g) = I dcrh{a)a'ui{qa/2)un{qcr/2). (A.13) 



The functions QiCn, ■ ■ ■ ,^n-i) are obtained from the same recurrence relations ([7]) with 
Ci = Poip^) defined now through the moments of the distribution functions. 
The system (lA.lOl) can be put in the matrix form 

e(o)(g) 
eW(g) 



M{po.q)\ ^(i)'J 1=0, (A.14) 



with 

Mio a)-I + o( '^'°'°'^^) ^^"''^"^ \ 6(^) - (e^'^ 6« Y (A 15) 

and / the 2n x 2n unitary matrix. The existence of a nontrivial solution of Eq. (1A.14P 
give us the following set of equations 

dM 
M{po,q) = -^{po,q) = 0, (A.16) 

(where we have defined M(po, q) = det [Ai{po, q)]) to find the values of Pq and q* (the 
absolute minimum of M{po, q) with respect to q) at bifurcation. 



Fluid-fluid versus fluid-solid demixing in mixtures of parallel hard hypercubes 18 

References 

[1] van Roij R, Mulder B and Dijkstra M, 1998 Physica A 261 374 

[2] Wensink H H, Vrocge G J and Lekkcrkcrkcr H N W, 2001 J. Chem. Phys. 115 7319 

[3] Dubois S and Pcrera A, 2002 J. Chem. Phys. 116 6354 

[4] Varga S, Galindo A and Jackson G, 2002 J. Chem. Phys. 117 7207 

[5] Schmidt M and Denton A R, 2002 Phys. Rev. E 65 021508 

[6] Martinez-Raton Y, Velasco E and Mederos L, 2005 Phys. Rev. E 72 031703; dc las Hcras D, 

Martincz-Raton Y and Velasco E, 2007 Phys. Rev. E 76 031704 
[7] Krancndonk W G T and Frenkel D, 1991 Mol. Phys. 72 679; 1989 J. Phys.: Condens. Matter 1 

7735 
[8] Dijkstra M, van Roij R and Evans R, 1998 Phys. Rev. Lett. 81 2268; 1998 Phys. Rev. E 59 5744 
[9] Almarza N G and Enciso E, 1998 Phys. Rev. E 59 4426 
[10] Martinez-Raton Y and Cuesta J A, 1999 J. Chem. Phys. Ill 317 
[11] Lafucntc L and Cuesta J A, 2002 Phys. Rev. Lett. 89 145701 

[12] Widom B and Rowlinson J S, 1970 J. Chem. Phys. 52 1670; Dickman R and Stell G, 1995 J. 

Chem. Phys. 102 8674; Brader J M and Vink R L C, 2007 J. Phys.: Condens. Matter 19 

036101; Schmidt M, 2004 J. Phys.: Condens. Matter 16 L351; Santos A, Lopez de Haro M and 

Yuste S B, 2010 J. Chem. Phys. 132 204506; Dijkstra M, 1998 Phys. Rev. E 58 7523 

[13] Conway J H and Sloane N J A, 1999 Sphere Packings, Lattices, and Groups, Springer, NY; Cohn 

H and Elkies N, 2003 Ann. Math. 157 689; Cohn H, 2002 Geom. Topol. 6 329 
[14 
[15 
[16 
[17 
[18 
[19 



[20' 
[21 
[22; 
[23; 
[24; 
[25; 

[26; 

[27 
[28; 
[29; 
[30 
[31 

[32; 
[33; 
[34 



van Meel J A, Frenkel D and Charbonneau P, 2009 Phys. Rev. E 79 030201 (R) 

Skogc M, Donev A, StilUngcr F H and Torquato S, 2006 Phys. Rev. E 74 041127 

Estrada C D and Robles M, 2010 arXiv:1002.3407v3 cond-mat.soft 

Finken R, Schmidt M and Lowen H, 2001 Phys. Rev. E 65 016108 

Robles M and Lopez de Haro M, 2007 J. Chem. Phys. 126 016101 

Bishop M and Whitlock P A, 2006 J. Stat. Phys. 126 299; Lue L, Bishop M and Whitlock P, 2010 

J. Chem. Phys. 132 104509 
Rohrmann R D and Santos A, 2007 Phys. Rev. E 76 051202 
Yuste S B, Santos A and Lopez de Haro M, 2000 Europhys. Lett. 52 158 
Lafuente L and Cuesta J A, 2002 J. Phys.: Condens. Matter 14 12079 

Cuesta J A and Martinez-Raton Y, 1997 Phys. Rev. Lett. 78 3681; 1997 J. Chem. Phys. 107 6379 
Tarazona P, 2000 Phys. Rev. Lett. 84 694 
Buhot A and Krauth W, 1999 Phys. Rev. E 59 2939 

Lafuente L and Cuesta J A, 2005 J. Phys. A: Math. Gen. 38 7461; 2004 Phys. Rev. Lett. 93 130603 
Tarazona P, Cuesta J A and Martinez-Raton Y, 2008 Lect. Notes Phys. 753 247 
Jagla E A, 1998 Phys. Rev. E 58 4701 
Groh B and Mulder B, 2001 J. Chem. Phys. 114 3653 

Hoover W G, Hoover C G and Bannerman M N, 2009 J. Stat. Phys. 136 715 
Hoover W G and de Rocco A G, 1962 J. Chem. Phys. 36 3141; Hoover W G and Poirier J C , 

1963 J. Chem. Phys. 38 327 
van Swol F and Woodcock L V, 1987 Mol. Simul. 1 95 
Kirkpatrick T R, 1986 J. Chem. Phys./ 85 3515 
Martinez-Raton Y and Cuesta J A, 1998 Phys. Rev. E 58 R4080 



